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Popular Summary 

Ideal cloud-resolving models contain so little accumulative errors that they can simnlatp 
cumulus clouds and synoptic large-scale circulations correctly. This paper discusses a discrete 
constraint for the motfels in relation to no accumulative errors or the correct representation of 
the large-scale convergence and advection of scalars (e.g. moist entropy, water vapor). The 
paper further shows that the challenge is how to compute sound waves efficiently under the 
constraint. 

To address the challenge, tire paper proposes a compensation method for the efficient 
computation of sound waves. Stability analysis and numerical experiments show that the 
method allows the models to integrate efficiently with a large time, step. 

In contrast to the time-split method, the compensation method needs no small-time-step 
integration for sound waves. As a result, the method is sli^tly more efficient than the time- 
split method when one processor is used, and much more efficient when many (e.g., hundreds 
of) processors are used for cloud-resolving modeling. 
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Abstract 


Ideal cloud-resolving models contain little’ acctunulative errors. When their domain 
is so large that synoptic large-scale circulations are accommodated, they can be used 
for the simulation of the interaction between convective clouds and the large-scale 
circulations, niis paper sets up a framework for the models, using moist entropy as a 
prognostic variable and employing conservative numerical schemes. The models 
possess no accumulative errors of thermodynamic variables when they comply with a 
discrete constraint on entropy conservation and sound waves. Alternatively speaking, 
the discrete constraint is related to the correct representation of the large-scale 
convergence and advection of moist entropy. 

Since air density is involved in entropy conservation and sound waves, the 
challenge is how to compute sound waves efficiently under the constraint. To address 
the challenge, a compensation method is introduced on the basis of a reference 
isothermal atmosphere whose governing equations are solved analytically. Stability 
analysis and numerical experiments show that the method allows the models to 
integrate efficiently with a large time step. 



1. Introduction 


Massive parallel computation upgrades the enviro nm ent for atmospheric modeling 
by domain decomposition (e.g., Droegemeier et al. 1995, Juang et al. 2003), breeding 
flie next generation of cloud-resolving models (e.g., a global cloud-resolving model). 
The future models will take so large a domain that they simulate both convective 
clouds and large-scale circulations explicitly, addressing the interaction among 
clouds, radiation and large-scale circulations, hi contrast to current cloud-resolving 
models, the future models will possess little accumulative errors for the correct 
simulation of large-scale circulations. 

Large-scale circulations are governed by thermodynamics in the Tropics 
(Raymond 1995, 2000; Raymond and Zeng 2000; Zeng, Tao and Simpson 2004). 
Zeng, Tao and Simpson (2004) used the analytical model of NeeUn and Held (1987) 
to show the sensitivity of tropical large-scale vertical circulations to atmospheric 
cooling rate. Their results indicated that, for the correct simulation of tropical large- 
scale vertical circulations, the accumulative temperature error should be much less 
than 1 K/day, the order of the atmospheric radiative cooling rate. In other words, the 
accumulative temperature error should be -10'^ K/day or less for the simulation of 
large-scale circulations. 

On the basis of the scale analysis of atmospheric convection (Ogura and Phillips 
1962), the current cloud-resolving models were constructed to simulate individual 
cloud systems for life-cycle characteristics (e.g., Klemp and Wilhelmson 1978, 
Grabowski 1989, Tao and Simpson 1993, Tompkins and Craig 1998, Xue et al. 2000, 
Tao et al. 2003). Some small terms were usually ignored for economical computation, 
which may distort the simulation of large-scale circulations. For example, the sink of 
moist air was ignored in expressing the air mass continuity equation in terms of the 



density of moist air. The approximation of no moist air sink can bring about an 
accumulative temperature error of 10'* K/day with a large-scale vertical velocity of 1 
cm/s (see Appendix A). Although such error can be ignored in the simulation of 
individual cloud systems for life-cycle characteristics, it can distort the simulation of 
tropical large-scale circulations. 

Moreover, numerical errors can distort the simulation of large-scale circulations. 
The current cloud-resolving models usually employ temperature (or its equivalent) as 
a prognostic thermodynamic variable. They compute sound waves economically with 
some approximations. If the approximations violate energy conservation, the 
accumulative temperature error may not be less than 10'* K/day. Consider a cloud- 
resolving model with a 10-second time step and a random temperature error around 
10'* K in one integration step. When the random error is averaged over many grid 
points and time levels, its value can exceed 10'^ K statistically in many cases, where 
the error of 10'^ K per step corresponds to an accumulative temperature error of 10'* 
K/day. Hence, large munerical errors can sometimes distort the simulation of large- 
scale circulations. 

The new trend is to replace temperature with moist entropy as a prognostic 
variable and diagnose the temperature from that and other prognostic variables 
(Raymond and Bl)dh 1986; Ooyama 1990, 2001; Zeng 2001; Zeng, Tao and Simpson 
2004). Moist entropy deals “easily” with the transition between the three water 
phases. If it is used as a prognostic variable, it can reduce numerical errors especially 
those connected with microphysical and dynamic processes (Zeng 2001; Zeng, Tao 
and Simpson 2004). Recently, Ooyama (2001) and Zeng (2001) constructed two- and 
three-dimensional models with moist entropy as a prognostic variable, respectively, to 
simulate warm clouds. Zeng, Tao and Simpson (2004) derived an accurate equation 
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for moist entropy, providing a theoretical basis for using moist entropy as a prognostic 
variable in long-term cloud modeling. Hence, one possible framework for the 
simulation of the large-scale circulation is a cloud-resolving model with moist entropy 
as a prognostic variable. 

If a cloud-resolving model uses moist entropy as a prognostic variable and 
discretizes the moist entropy and other prognostic thermodynamic variables in 
conservative forms, it contains no accumulative error of the prognostic 
thermodynamic variables. Since temperature is diagnosed from the moist entropy and 
the other prognostic thermodynamic variables, the model contains no accumulative 
temperature error either. However, the conservative discrete form of the moist entropy 
involves air density, and the air density is related to sound waves through the air mass 
continuity equation. Although sound waves are meteorologically unimportant, their 
phase speed limits the time step for integration. Thus, the challenge is how to compute 
sound waves economically while conserving moist entropy (see Section 2 for more 
discussion). 


Many good n umerical schemes for sound wa v es h a ve b e en proposed for life - cy ele- 
modeling of individual cloud systems (e.g.. Hemp and Wilhehnson 1978, Anderson, 
et al. 1985, Droegemeier and Wilhehnson 1987, Skamarock and Hemp 1992) and 
long-term cloud ensemble modeling (Zeng 2001, Satoh 2002). However, those 
schemes do not suit the cloud-resolving model well for the simulation of large-scale 
circulations, because the model is expected to possess the following three 
characteristics: (1) efficient computation with a single processor and highly efficient 
massive parallel computation with 10 -10 processors, (2) entropy or energy 
conservation, and (3) no accumulative error even from the computation of sound 


waves. To prepare such models for the simulation of large-scale circulations, this 
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paper introduces an efficient method for computing sound waves on the basis of a 
reference atmosphere whose governing equations are solved analytically. 

With the explicit simulation of convective clouds and large-scale circulations as 
the background, this paper proposes a framework for future cloud-resolving models 
and deals with the computation of sound waves in the finmework. The paper consists 
of five sections. Section 2 introduces a framework for cloud-resolving models with 
moist entropy as a prognostic variable. It explains a discrete constraint on entropy 
conservation and sound waves that corresponds to no accumulative error of 
thermodynamic variables. Section 3 develops an efficient method to compute sound 
waves, and Section 4 tests the method with numerical experiments. Section 5 gives a 
summary. 

2. Entropy conservation and sound waves 

Accumulative errors in a numerical model originate in the continuous governing 
equations and their discretization. When the air mass continuity equation is expressed 
in the density of moist air and the sink of moist air due to water vapor condensation is 
ignored, there is an accumulative temperature error in the continuous governing 
equations. The error, as shown in Appendix A, is so large that it sometimes can distort 
the simulation of tropical, large-scale circulations. To remove the error and similar 
others, a finmework for cloud-resolving models is set up for the simulation of large- 
scale circulations. 

a. A framework for cloud-resolving models 

A framework is set up for cloud-resolving models with no accumulative error of 
thermodynamic variables, where the density of dry air and the moist entropy are used 
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as prognostic variables. Specifically speaking, the density of dry air /> is used to 
rewrite the ideal gas law, the air mass continuity equation and the flux form of 
prognostic thermodynamic variables. 

Ibe ideal gas law for moist air is written as 

p = pRJ(X + l.m%q^) (2.1) 

where p is the total pressure of moist air, T the air temperature, qy the mi xin g ratio of 
water vapor and Rd the gas constant of dry air. 

Tbe momentum equation for moist air (including precipitating particles) is written 
as 

d,v + /jqxv + VAf = — Vp + g + D, (2.2) 

Pil + q.+q^) 

where v is the velocity vector and g the acceleration due to gravity; Dv represents the 
subgrid diffusion of momentum due to turbulence; tT = -j v • V is the kinetic energy 
per umt mass; q = yO ’(Vx v + 2£l) is the vorticity divided by air density where £2 is 
the angular velocity vector of the Earth’s rotation; q, is the total mixing ratio of 
airborne water and qp the mixing ratio of precipitating particles. If qc, qt, qr, qs and qg 
denote the mixing ratios of cloud water, cloud ice, rainwater, snow and grauple/hail, 
respectively, the total mixing ratios of airborne water and precipitating particles are 
expressed as 

9r=9v+9c+9i 

The mass continuity equation for diy air is 

a,yO + V-(yOv) = 0 (2.3) 

and the prognostic equation for a scalar (p \s 
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h,p<t> + V ■ {ps$) = p{M^ + D^) 


(2.4) 


where refers to the microphysical processes in clouds and the scalar subgrid 
mixing due to turbulent eddies. 

The sjiTnbol Equation (2.4) represents s the moist entropy, qu qu qr, qs and qg 
the mixing ratios of airborne water, cloud ice, rainwater, snow and graupel/hail. 
Following the notation of Tao and Simpson (1993), the microphysical processes are 
described by q, , qi, q^, qs and qg) or 


M^ = £, + S,-t5,-(T,,+r„.) 

(2.5a) 

•9 

II 

1 

1 

9^ 

(2.5b) 

M^=p-'a,(pV,q,)-E,-F,-F, -r„ 

(2.5c) 


(2.5d) 


(2.5e) 

M, = a, + (Q-yD,)/T 

(2.5f) 


where F, is the terminal velocity of precipitating particles; E, F and 5 stand for 
evaporation, fusion/freezing and sublimation/deposition, respectively; and Tgc, Tgr, Tgi, 

and ^98 represent the microphysical transfer rates between hydrometeor species 
and their sum is zero. 

In Equation (2.5f), Q is the rate of diabatic heating and Gs the internal source of 
moist entropy; the term -v-D^ represents the heat generated by internal friction and is 
introduced for energy conservation. The moist entropy s is defined iu unit mass of dry 
air, representing the contribution from dry air, water vapor, cloud water and ice except 
for precipitating particles. The complete flux-form equation for moist entropy was 
introduced by Zeng, Tao and Simpson (2(X)4). In other words, the expression for Gs is 
known. 
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The models employ the following prognostic variables: the velocity vector v, the 
air density p, the moist entropy s, the total mixing ratio of airborne water qt, the 
mixin g ratio of cloud ice qi, and the mixing ratios of rainwater qr, snow qs and 
graupel/hail qg. Their corresponding prognostic equations are Equations (2.2), (2.3) 
and (2.4). Once the prognostic variables are known, the total pressure of moist air p is 
determined by Equation (2.1), and the quantities T, qv and qc are diagnosed from the 
prognostic variables s, qt and qi. 

b. A discrete constraint 

Scalars are conserved in advection, to which the intension of “conservation” is 
confined in this paper. It is easy to apply conservative schemes to Equations (2.3) and 
(2.4). When the scalar (pis constant and the discrete form of Equation (2.4) 

should reduce into that of Equation (2.3), which is referred to here as a discrete 
constraint. Recently the constraint has attracted attention because of its importance in 
the construction of a model for the interaction between tropical convection and large- 
scale circulations (Zeng 2001, Arakawa 2004). Zeng (2001) constructed a three- 
dimensional cloud-resolving model under the constraint to simulate tropical 
convective clouds and their interaction with large-scale circulations. Arakawa (2004) 
referred to the constraint as “constancy” warranty and related the constraint to the 
correct representation of the advection and convergence of scalars such as water 
vapor. 

When the models in the preceding subsection employ conservative schemes, the 
constraint corresponds to no accumulative error of thermodynamic variables. In the 
models without the constraint, accumulative errors exist, one of which on p is 
understood from the following mathematical operations. Substimting ^in the discrete 
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form of Equation (2.4) with^4:» then subtracting the discrete form of Equation (2.3) 
times <j>c yields a resulting equation, where the constant ^ is chosen arbitrarily. In 
comparison witii the discrete form of Equation (2.4), the resulting equation contains a 
new string that is composed of all terms with (4- The new string represents an error 
and is proportional to ^c- When the error is accumulated, the simulation of large-scale 
circulations may be distorted^ 

In fact, it is easy to discretize Equations (2.3) and (2.4) under the constraint. A key 
question is how to compute sound waves efficiently under the constraint. Since the 
discrete form of Equation (2.3) is related to sound waves through the ideal gas law, 
the speed of sound waves limits the time step for integration. To construct an efficient 
scheme for sound waves, the ideal gas law or Equation (2.1) is modified to 

1-608 p-A?-R^r(l + 1.608qJ 

Bt ^^pdt Tdt 1 + 1.608^, dt^ T, ^ ^ ^ 

where the relaxation timescale zj is chosen to be 5 minutes or so. 

Equation (2.6) is close to the ideal gas law. Once it deviates from the ideal gas law 
as it sometimes can due to computational reasons, it soon relaxes to the ideal gas law 
on the timescale Its temporary deviation from the ideal gas law provides room for 


* For a given numerical scheme, the error is estimated as follows. Assume that ^ is constant. 
Dividing the discrete form of Equation (2.4) by then subtracting the discrete form of 
Equation (2.3), and finally moving all terms to the right side yields a resulting equation. The 
right side equals zero under the constraint. Otherwise, the right side represents a spurious 
mass sink in a discrete equation for either air mass continuity or scalar advection. The 
spurious mass sink is compared with the sink of moist air due to water vapor condensation, 
estimating the order of the error just as doing in Appendix A. 
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the numerical treatment of sound waves under the constraint and brings about no 
accumulative error of thermod 3 mamic variables. 


3. A compensation method for sound waves 

An efficient method is introduced in this section to compute sound waves under 
the constraint. For srmplicity, the mediod is illustrated in a one-dimensional dry 
model as an example without loss of generality. From Equations (2.2), (2.3) and (2.6), 
the equations for a one-dimensional dry atmosphere are written as 


d,pw + d.p + pg+- = 0 

(3.1a) 

d,p + d^pw + --- = 0 

(3.1b) 

d,p- RJd,p + ip - pRJ)r;' + — = 0 

(3.1c) 

where z is the height and w the vertical velocity component; and the temperature T is 

fixed for the convenience of description. 


In contrast to Equation (3.1), a simple reference isothermal 

atmosphere is 

introduced with the following governing equations 


d,pw + d^p + pg =0 

(3.2a) 

djP + d ^pw = 0 

(3.2b) 


(3.2c) 


where the temperature Tref is a constant Subtracting the left sides of Equation (3.2) 
from both sides of Equation (3.1) yields 


{d,pw + d^p + pg+---) 

+ Pg)=-^,{pwY‘‘^ +pg) 


(3.3a) 


{d,p + d^pw+--) 

-(a, +a,(pw)^‘'>)=-(a,p<“> +a^(/Av)^“>) 


(3.3b) 
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(3.3c) 


(a,p - RJd,p + (p - pRJ)tf +■■■) 

where the superscripts {d) and (a) are introduced only to group terms for different 
numerical schemes. Strictly speaking. Equation (3.3) is equivalent to Equation (3.1) 
because a variable with the supCTScript (^^) represents the same property of die real 
atmosphere as its corresponding variable with the superscript (a). Hence 


d,pw = a, (/xv)'*'^ = a, (/xv)^“> 

(3.4a) 

^,p = d,p^‘^'> ^d,p^“^ 

(3.4b) 

a,p = ay‘'» = a,p<“> 

(3.4c) 


a. Symbol notation 

Equations (3.1), (3.2) and (3.3) are discretized on three time levels, from the time 
level n-1 to n+1, and four sets of symbols are introduced for concise expressions. 
Any variable for a given property in the three equations takes the same values on the 
time levels n-1 and n except for on the time level n+1. The symbols 
p(n+D) denote the final values of variables on the time level n+1 in the numerical 
method, and (w*, /?*, p*) the values on the time level n+1 that are calculated from the 
finite-difference form of Equation (3.1). If (w*, /?*, p*) is set, the 

time step for integration is limited by die speed of sound waves, which degenerates 
into a traditional method. 

The symbols (w^*®, yC^‘®, p^"^) and (>v^"\ denote the two sets of variables on 

the time level n+1. They are calculated numerically from the finite-difference form of 
Equation (3.2) and analytically from Equation (3.2), respectively. Because (w^'^, 
p^“^) and (w^“\ ff°\ p^“^) approximate the same properties of the reference atmosphere 
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2 

in two different ways , their difference is close to zero theoreticaUy. Thus the values 

still approximate the solution of Equation (3.1) 
or (3.3) on the time level n+1. In other words, con^nsate 

(w*, p, p*) for the final values, which increases the maximum time step for stable 
integration and is referred to here as tibe compensation method. 

The four sets of symbols are illustrated in the preceding paragraphs. Next is the 
computational procedure in the compensation method with accurate definitions of the 
symbols. Apply similar finite-difference schemes to Equations (3. 1) and (3.2), and use 
their discrete forms to calculate (w*, p*, p*) and (w^*^, p^^ on the time level «+l, 

respectively. Then introduce the symbols 


Sc.=((Av)'"^‘-(pw)<"'"*‘)/2A/ 

(3.5a) 


(3.5b) 

S^=(p‘''"’-p'‘'’"*')/2Ar 

(3.5c) 

where At is the time step and the superscript n+1 is added to distinguish the variables 

on the time level n+1 from the variables on other time levels. Using Equations (3.4) 
and (3.5), Equation (3.3) is rewritten as 

d,(pw)^“>+a,p'‘'^+/7'-g = 5^ 

(3.6a) 

d,p^“Ud^(pwr^=s,^ 

(3.6b) 


(3.6c) 


^ The symbols (w^'^, p*-^ and (w^“', refer to the reference atmosphere only in this 

paragraph for easy understanding. They are defined accurately in the following paragraphs. In 
fact, it is unnecessary for a mathematician to introduce the reference atmosphere because the 
method does not need an equation like Equation (3.2). 
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with the initial condition 




= p"-> 


(3.6e) 
(3.6f) 

at t={n-l)At and die boundary condition 

= p, atz=o (3.6g) 

(pH')^“^=0 at z=z, op, (3.6h) 

where Zu>p is the height of the top boundary, and the variables with the subscript “1” 
refer to the properties near the surface. 

Equation (3.6) is solved analytically for the values at 

r=(ra+l)Af. From this and Equation (3.4) the final solution of Equation (3.1) at the 
time level n+1 is obtained. That is 

(3.7b) 


b. Solution expression 

This subsection shows how to solve Equation (3.6) analytically. Let Cs,={RdTref)^'^ 
denote the sound wave speed in the reference atmosphere and Zs=lc„tsi the distance 
that the sound waves travel in 2Ar. Assume that Sew, Sep and Sep do not change with 
time between the time levels n-1 and n+1. From the boundary condition of Equation 
(3.6), w and p are set as 


-ip (z) + S^^ (z))dz + (c;^5^(z) + 5,^ (z))dz 


(3.8a) 
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+ p:-S |>”* 

.fz+z,/ , V ° „ (3-8b) 

+ 2‘^«-J (■^q,(z) + c^S (z))rfz-|f ‘S^{z)dz + \\ S^{z)dz 

•*i i, JZ Jt-Zf 


where ps is the surface pressure; and the new variables W\z, t) and p\z, t) satisfy 

a,W(2,0 + 9,p'fe0=0 (3.9a) 

a,p'(z,t) + c^a,W'(z,O = 0 (3.9b) 

wife fee boundary conditions W (z, t)=0 at z=0 and Zwp, and fee initial conditions 

W\z,t) = W^(z) = p"-^w'-^ -pr'<'(X-z/z,^) (3.9c) 


p'iz,t) = Poiz) = p"-' -p"+g 
at /=(w-l)Ar. Equation (3.9) is solved analytically, giving 

+ Po(z + zJ- Po(z-Z,) 

2 2c^ 

. + + Wo(z + zJ-Wq(z-z,) 

^ ^ ^sr ^ 


(3.9d) 


(3.10a) 

(3.10b) 


at t=(n+l)Ar, where the function Wo and po are extended to (-oo,+oo) wife their 
symmetry to z=0 and Ztop- 

To calculate w and p on fee time level n+1 wife Equations (3.8) and (3.10), fee 
values of w and p on fee time level n-1 between grid points are needed. They are 
calculated by linear interpolation. As a side effect of fee linear interpolation, sound 
waves are damped. 

The air density is calculated wife Equation (3.6c) after the pressure p"*^ is 
known. Since fee pressure is calculated with the aid of linear interpolation, fee 
vertical integration of ff*^-p may not equal zero. For mass conservation, a slight 
adjustment oi is taken so feat the vertical integration of equals zero. 
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c. Stability explanation 

The numerical stability of the compensation method is analyzed in Appendix B for 
a linear system that is close to the preceding one. Here is a physical explanation for 
the numerical stability. The method, as Equations (3.3) and (3.6) show, works like 
two systems. The two systems step forward in turn. One system is governed by 
Equation (3.3) with zeros on its right sides. Its sound wave speed (Rj(T is 

low. It is simulated with finite-difference schemes. Its time step is limited by the low 
sound wave speed if no other system is involved. 

The other system is governed by Equation (3.6). Its sound wave speed 

is high. The system is simulated with anal)rical schemes instead of finite-difference 
schemes. Thus the time step for integration is not limited by the high sound wave 
speed. 

The two systems step forward in turn. First the slow “sound waves” propagate 
locally. Then their residue is carried by the fast “sound waves” into a vast space. As a 
result, the “sound waves” accomplish the adjustment between the velocity and the 
pressure fields in the models even though the “sound waves” do not imitate real sound 
waves well. Since the time step for integration is not limited by the high sound wave 
speed, the models can be integrated with a large time step. The maximum time step 
for stable integration can be determined by the Fourier stability analysis such as in 
Appendix B or the numerical tests in the next section. 

4. Numerical tests 

To test the compensation method, a one-dimensional atmosphere that is governed 
by Equation (3.1) without ellipses is chosen as a test case. Initially the atmosphere is 
stationary and in hydrostatic equihbrium; the surface pressure is 1013.25 hpa; the 
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temperature decreases linearly with height from the surface value of 288 to 216.5 K at 
z=ll km and remains constant above 11 km. The atmospheric rigid boundaries are 
located at z=0 and Zjop=18 km. A vertical force 

F, = 10~^ g sin(— ) sin (— — ) when z<2ztoJ'i 

is applied to a unit mass of air, where the timescale minutes. Hie force, 

imitating the buoyancy in convective cells with updrafts and downdrafts in turn, 
drives the atmospheric motion. 

Two experiments with a vertical grid size Az=600 m are done to simulate the 
atmosphere, using the compensation method and the traditional method in Appendix 
B, respectively. A time-smoother (Robert 1966) is used to remove any tendencies that 
might decouple the odd and even steps. In the experiment with the traditional method, 
the time step A/=0.1 seconds altiiough the maximum time step for stable integration is 
close to 1 second. Since the time step is so small, the results in the experiment are 
treated as a standard benchmark to measure the results in the other experiment. One 
result in the experiment, the vertical velocity pwlps at z=6 km with Pj=l kg/m^ versus 
time, is displayed in the upper panel of Figure 1, showing the superimposition of the 
oscillations of sound waves on a forced oscillation. 

In the experiment with the compensation method, tiie time step An=10 seconds, the 
reference temperature Tr^273.15 K, and the relaxation timescale ts=5 minutes. The 
vertical velocity at z=6 km versus time is displayed in the lower panel of Figure 1. Ih 
contrast to the results in the experiment with the traditional method, the forced motion 
due to buoyancy is simulated well although the sound waves are not simulated well. 

The two experiments are re-done to show why the sound waves are not simulated 
well in the compensation method. In the two new experiments, a sponge layer is 
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introduced above 15 km to absorb sound waves. The Tnaximiim time scale for 
damping sound waves is 30 seconds at the top boundary. Figure 2 displays the same 
variable in Figure 1 for the new experiments. As the figure shows, the forced 
oscillations in the two experiments are similar. In the traditional method, the force 
generates sound waves and the sound waves are damped in the sponge layer. In 
contrast, in the compensation method, not only the force but also the numerical 
schemes generate sound waves. The sound waves are damped not only in the sponge 
layer but also by the numerical schemes. Based on Figure 1, the sound waves are 
mainly damped by the numerical schemes. Otherwise, the amplitude of sound wave 
fluctuation would increase with time. 

Since the numerical schemes in the compensation method generate sound waves 
continuously as computational noise, an interesting question is whether the noise 
affects the simulation of meteorological phenomena. If the numerical scheme for 
scalar advection in Appendix C is used corresponding with the compensation method, 
sound waves bring about no accumulative error of thermod 5 mamic variables. As a 
result, sound waves do not affect meteorological motion. 

Two experiments are conducted to show that sound waves do not affect 
meteorological motion. In one of them, the traditional method with a small time step 
(Ar=0.1 s) is used, and a sponge layer for sound waves is introduced. Results for the 
first five hours are shown in the upper panel of Figure 2, and the results after 70 days 
are shown by the thick line in Figure 3. Since sound waves are damped completely in 
70 days, the thick line represents the meteorological motion. In the other experiment, 
the compensation method with a large time step (Ar=10 s) is used, and there is no 
sponge layer for sound waves. Results for the first five homs are shown in the lower 
panel of Figure 1, and the results after 70 days are shown by the thin line in Figure 3. 
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Figure 3 clearly shows that the compensation method simulates the meteorological 
motion quite well in the long-term integration in the presence of computational sound 


waves. 

5. Discussion and summary 

On the basis of the scale analysis of atmospheric convection (Ogura and Phillips 
1962), almost all current cloud-resolving models ignored some small terms for 
economical computation. Those models have successfully simulated convective cloud 
systems for life-cycle characteristics. When they are extended to the simulation of 
large-scale circulations, however, the question is whether the small terms are still 
ignorable. 

Zeng, Tao and Simpson (2004) used the analytical model of Neelio and Held 
(1987) to diagnose large-scale vertical velocity in the Tropics. Their results showed 
the sensitivity of large-scale vertical circulations to the atmospheric coohng rate and 
the surface flux of moist entropy from the underlying surface to the air above. Since 
the atmospheric radiative cooling rate is 1 K/day or so, their results indicate that the 
accumulative temperature error should be equal to lO ’ K/day or less for the 
simulation of tropical large-scale vertical circulations. 

Using the value of 10 * K/day as a scale, the approximations in the current cloud- 
resolving models are re-checked for the simulation of large-scale circulations. It is 
found that some approximations should be removed. For example, the density of 
moist air should be replaced with the density of dry air in expressing the air mass 
continuity equation (see Appendix A), and the thermal capacity of precipitating 
particles should be taken into account (Zeng, Tao and Simpson 2004). A discrete 
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constraint on the consistency between the discrete scalar equation and the discrete air 
mass continuity equation should be addressed. 

For the correct simulation of large-scale circulations, a framework is set up for 
cloud-resolving modds that use the density of dry air to express the air mass 
continuity equation and employ moist entropy as a prognostic variable. In the models, 
the flux-form equation for moist entropy takes account of the thermal capacity of 
precipitating particles (Zeng, Tao and Simpson 2004). 

The models possess no accumulative error of thermodynamic variables when they 
comply with a discrete constraint on entropy conservation and sound waves. The 
constraint requires that the discrete flux -form equation for scalar advection reduce 
into the discrete equation for air mass continuity while the scalar is constant. Since air 
density is involved in the flux-form equation for scalar advection and is related to 
sound waves through the air mass continuity equation, the challenge is how to 
compute sound waves efficiently under the constraint. 

To address the challenge, a compensation method is developed for the efficient 
computation of sound waves, where a reference isothermal atmosphere is introduced. 
The governing equations for the reference atmosphere are solved both analytically 
and numerically. The difference between the analytical and the numerical solutions is 
used to compensate the original solution of the models with traditional finite- 
difference schemes. Once the original solution is compensated, the models simulate 
sound waves equivalently by “analytical” schemes rather than finite-difference 
schemes (see the text for an accurate description). As a result, the time step for 
integration is not limited by the high speed of sound waves, and the models can 
integrate efficiently with a large time step. 



The compensation method for soimd waves is tested in a one-dimensional model. 
Numerical experiments show that the method with a large time step is stable and can 
simulate meteorological motions quite well in a long-term integration. On the same 
principle, the method can be extended to multi-dimensional models easily. The 
compensation method in a three-dimensional, terrain-following coordinate system is 
technically complicated and will be introduced elsewhere as a technical report. 

In contrast to the time-split method (Klemp and WiUiehnson 1978), the 
compensation method needs no small-time-step integration for sound waves. Thus, 
when one processor is used for computation, the models with the compensation 
method ran slightly faster than those with the time-split method. Recent parallel 
computation tests of the Goddard Cumulus Ensemble model (Tao et al. 2003) on 
memory-distributed machines showed that the efficiency of parallel computation is 
degraded by the data communication between processors in computing sound waves, 
especially while hundreds of processors are used. Without smaU-time-step integration 
for sound waves, the compensation method is expected to be much more efficient than 
the time-split method in massive parallel computation for cloud-resolving modeling. 
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APPENDIX A 


Analysis of Accumulative Temperature Error 
In this appendix, the temperature error is analyzed in a set of continuous 
governing equations, one of which,.the air mass continuity equation, ignores the sink 
of moist air due to water vapor condensation or deposition. To show the error, the 
density of moist air /C% is distinguished from the density of dry air pd by 
/>, = (1 + ^„ ) . Differentiating the proceeding relation yields 

dlapj dt ~ d\a p^l dt + dq^l dt (Al) 

If the potential temperature 0 is described by dOjdt = Mg where Me represents the 


source of potential temperature, the accurate flux form of <9 is 

dp.d 

J^^^.{p^6v) = p,Mg (A2) 

When a cloud-resolving model uses the following prognostic equation 

+ V • ip^Ov) = p,Mg (A3) 

it possesses a temperature enor. The mathematical operation, multiplying Equation 
(A2) by pjpd and then subtracting Equation (A3), yields the expression for the error 

^ dt , dt 

with the aid of Equation (Al). 

In cloudy air, qv=qvs the saturation mixing ratio of water vapor. Using the 
Clausius-Clapeyron equation and the hydrostatic balance approximation, the 
proceeding expression for the error is simplified into 




(A4) 
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where w is the vertical velocity, ys the temperature lapse rate of saturated air, Ly the 
latent heat of vaporization and the gas constant of water vapor. Let 0~T -273 K, 
qy,~5 g/kg, and >5-0.6x10'^ K/m. Thus ~ 10'’ K/day when w=0.01 m/s. 

The error increases with the vertical velocity. When w=0.1 m/s, the error is 10° K/day, 
which is the order of the atmospheric radiative cooling rate. 

Zeng, Tao and Simpson (2004) used the model of Neelin and Held (1987) to show 
the sensitivity of the large-scale vertical circulations to the atmospheric cooling rate in 
the Tropics. Their results indicated diat the accumulative temperature error should be 
~ 10 * K/day or less for the correct simulation of tropical large-scale vertical 
circulations. Hence the approximation in Equation (A3) should be avoided in the 
numerical simulation of tropical, large-scale circulations. 
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APPENDIX B 


Fourier Stability Analysis 

The computational stability of the compensation method is analyzed in this 
appendix. Similar to the system in Section 3, the following linear system is introduced 


= 0 

(Bla) 

0,/> + 0^W = O 

(Bib) 

^tP-c]'^tP^ip-c]p)lT:, =0 

(Blc) 


where (w, p, p) are prognostic variables, and (c^, T,) are constants. It is discretized on 
three time levels as 


^ *«+l .,n-l f n \ 

^ (Pj+l/2 ~Py-l/2) 

^n+1 yn"”! 2Ar , I, „ . 

Pj+V2 ~ Pj+V2 ~ ^ 


(B2a) 

(B2b) 


= p;:L - p;:u2)+—(p-:12 -c^pj:L) ob2c) 

where At is the time step and Az the spatial grid size; the variables w and p are 
spatially staggered; and the variables with an asterisk will be used to calculate their 
final values on the time level n+1. If the variables with an asterisk are treated as their 
final values. Equation (B2) is a traditional method where the time step is limited by Cs 
the sound wave speed. 

To increase the time step for integration, the following reference system is 
introduced just as in Section 3 


(B3a) 

(B3b) 


d,w + d^p = 0 

d,p + dpv = Q 

^tP-cld,P = 0 
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(B3c) 


where c„ is the sound wave speed in the reference system. Equation (B3) uses the 
same finite-difference scheme as Equation (B2). After the different values of the 
variables on the time level n+1 in the two systems are obtained, the expressions for 
5cw and Sep in Equation (3.5) become 

(B4a) 

Wu 2 ) = -ic] -cl)(w%, -w])/Az + -c]p%^)/T^ (B4b) 

From Equations (3.8) and (3.10), the final values of the variables on the time level 
w+1 are obtained. For the sake of convenience, it is assumed here that rj’ = 0 and 


•^cp(j+l/2) ~ ^ 


M — m- 


2c, 


(B5) 


where m=l, 2, ... is an integer. With the use of Equations (B4) and (B5), w"""’ and/?"'^^ 
are expressed firom Equations (3.8) and (3.10) as 


,B+I _ ^i-m _ (Pjtl/2-Kn Pj-V2*m ) (P j+V2-m P j-J/2-m ) 

2 4c„ 


,.n-l 


W 


(B6a) 


„n+I _ Pj+U2+m + Pj+V2-m + ^"+1 ) “ ) 

Pj-HI2 ~ ^ 

2 4 

‘■n-v 2 A 

<^sr 4 


(B6b) 


Fourier stability analysis (e.g., Lomax et al. 2001) is applied to Equation (B6). 
Consider a harmonic exp(c» + ikz) where k is real. Let a = exp(oAr) . Thus |c^<l for 
numerical stability. Substitution of the harmonic into Equation (B6) yields 

(cr^ -cos(kmAz))^ + a{c]c^ -l)(l-cos(fcmAz))(( 7 ^ -cos(kTnAz)) 


+ cos (Mz / 2) sin ^ (kmAz)(l - cr(c^cj - 1)) = 0 
When Cs,=Cs, the expression for cris obtained from Equation (B7), or 


(B7) 
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\(t\ —l-sin^(kAz/2)sm^(kmAz) (gg) 

Thus |o|<l which shows that the analytical solution is stable. All modes decay in time 
except for the modes with sin(ibnAz) = 0 . 

When and sin(fcmAz) = 0, the expression for cris obtained from Equation 
(B7), or 


1 ct |=1 


(B9) 


when 1 1 [< 1 . Equation (B9) shows that the modes with sin(;tmAz) = 0 do not 

grow in time when 
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APPENDIX C 


A Numerical Scheme for Scalar Advection in the Compensation Method 

A numerical scheme is introduced in this appendix for scalar advection under the 
discrete constraint. First the discrete air mass continuity equation in the compensation 
method is summarized. Then a scheme for scalar advection is proposed on the basis of 
the discrete air mass continuity equation. 

In the compensation method, there is a discrete form of the air mass continuity 
equation 

/^j+l/2 ~Pj+l/2 ^ ^Pi+^n'^ Pj+V2)^j+l~^Pj^\n'^ Pj-UTy^J ^ ... _Q (Cl) 

2 Ar 2 Az 

where Az is the spatial grid size and the variables w and p are spatially staggered; the 
superscript n and the subscript j indicate time level and space grid-point, respectively; 
the air density p*"^2 > a temporary value, is different from the final value ™ 
compensation method. 

Since the time step Ar is large in the compensation method, the distance that sound 
waves travel in Ar is much larger than Az. Thus the cloud-resolving model is 
computationally unstable if p^j/2 = P/+1/2 is in other words, is determined 


by not only and but also the vertical velocity at other grid-points. 

In the compensation method, p"+]/2 is known at last. With reference to Equation 
(Cl), ppj,2 is expressed using w" and in the following form 

2 Ar 2 Az 

(C 2 ) 
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where the equivalent velocity is introduced to explain ~ *e 


deviation of the final air density from the prediction of Equation (Cl) at 2At due to 
sound waves. Subtracting Equation (Cl) from (C2) yields 

2Ar 2Az “ 

which determines with the aid of the boundary condition at z=0 and Ztop- 

Obviously w^^^is much smaller than the phase speed of sound waves. 

Once is known, it is easy to construct a discrete form of the equation for 
scalar advection under the constraint. With reference to Equation (C2), Equation (2.4) 
can take the following discrete form 


Hj^V2Vj+V2 


2Ar 


- + 


(P%3n<P%3n+P%U2Ku2)(^%. +^%l-iP%u2<^Un+Plu2<^lv2)i<+wf’‘) 


2Az 


= 0 
(C4) 


When the scalar <^is constant. Equation (C4) reduces into Equation (C2). In other 
words, the scheme in Equations (C2) and (C4) complies with the discrete constraint, 
and Equation (C4) does not change scalar advection. 
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Figure Captions 


Figure 1 The vertical velocity at z=6 km changes with time in the two experiments 
with no sponge layer for sound waves. Hie results in the traditional and the 
compensation methods are shoivn in the upper and the lower panels, respectively. 

Figure 2 the same as Figure 1 except that a sponge layer for sound waves is 
introduced. 


Figure 3 The vertical velocity at z=6 km changes with time in two experiments. The 
thick and the thin lines show the results in the two experiments with the traditional 
and the compensation methods, respectively. In the experiment with the traditional 
method, A/=0.01 s, and a sponge layer for sound waves is introduced. In the 
experiment with the compensation method, At=10 s, and no sponge layer for sound 
waves is used. 
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Figure 1 The vertical velocity at z=6 km changes with time in the two experiments with 
no sponge layer for sound waves. The results in the traditional and the compensation 
methods are shown in the upper and the lower panels, respectively. 
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Hjgure 3 TTie vertical velocity at z=6 km changes with time in two experiments. The 
thick and the thin lines show the results in the two experiments with the traditional and 
the compensation methods, respectively. In the experiment with the traditional method, 
A^=0.01 s, and a sponge layer for sound waves is introduced. In the experiment with the 
compensation method, At=10 s, and no sponge layer for sound waves is used. 
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